The Interplay of Nonlinearity and Architecture in Equilibrium Cytoskeletal Mechanics 
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The interplay between cytoskeletal architecture and the nonlinearity of the interactions due to 
bucklable filaments plays a key role in modulating the cell's mechanical stability and affecting its 
structural rearrangements. We study a model of cytoskeletal structure treating it as an amorphous 
network of hard centers rigidly cross-linked by nonlinear elastic strings, neglecting the effects of 
motorization. Using simulations along with a self-consistent phonon method, we show that this 
minimal model exhibits diverse thermodynamically stable mechanical phases that depend on ex- 
cluded volume, crosslink concentration, filament length and stiffness. Within the framework set by 
the free energy functional formulation and making use of the random first order transition theory 
of structural glasses, we further estimate the characteristic densities for a kinetic glass transition 
to occur in this model system. Network connectivity strongly modulates the transition boundaries 
between various equilibrium phases, as well as the kinetic glass transition density. 



I. INTRODUCTION 

The cytoskeleton is a crowded network of dynamic fil- 
amentous protein polymers that collaborate with diverse 
binding proteins and molecular motors to form nature's 
most marvelous active material^. All eukaryotic cells are 
known to contain a well-developed cytoskeleton, and even 
bacteria have been found to contain a diverse set of pro- 
teins capable of forming structural filaments 2 -. Actin is a 
major constituent of the cytoskeletal network that partic- 
ipates in such widely differing processes as blood clotting, 
egg fertilization, intestinal absorption and tumor inva- 
sion. A common theme of many substantial experimen- 
tal efforts^ — is the role of actin filaments in maintaining 
cell architecture and generating movement. 

In vitro, actin can polymerize to form long rigid fila- 
ments (F-actin), with a diameter around 7 nm and con- 
tour length up to 20 fim. The in vivo cytoskeletal net- 
work, however, is regulated and controlled not only by 
the concentration of F-actin, but also by accessory pro- 
teins that bind to F-actin. Nature provides a host of 
actin-binding proteins (ABPs) with versatile functions 
that offer the necessary variations for the part actin 
has to play^: Cross-linking proteins form filament bun- 
dles or isotropic gels, whereas capping proteins and/or 
filament-severing proteins regulate actin polymerization 
under specific salt conditions and thus control the average 
length of F-actin. Experiments show, for example, that 
actin-binding protein functions to stabilize cortical actin 
in vivo and is required for efficient cell locomotion^. Rhe- 
ological experiments on simplified networks and living 
cells demonstrate that even the most basic mechanical 
properties of cytoskeletal materials are sensitive to spe- 
cific architectural details, including filament and cross- 
link density, connectivity, and orientation^. 

Unlike flexible polymers, where changes in cross-link 
density typically do not markedly affect the elastic- 
ity, small changes in the concentration of cross-linking 



actin-binding proteins do dramatically alter the elastic- 
ity of F-actin networks^. As materials, in vitro net- 
works of cytoskeletal filaments exhibit several unusual 
mechanical properties, including a highly nonlinear elas- 
tic respons o 7 ' 12 and negative normal stresses^. Another 
noteworthy observation is the buckling of actin stress 
fibers that occurs upon rapid shortening^, which demon- 
strates the instability of a prestressed actin network un- 
der compression. 

A refined experimental route to understanding cy- 
toskeletal structures in terms of their molecular compo- 
nents is to reconstitute these structures from purified pro- 
teins. In vitro reconstitution provides the capability to 
study the emergence of micrometer-scale function from 
numerous molecular-scale interactions. A recent recon- 
stitution of contractility in a simplified model system, 
composed of purified F-actin, muscle myosin II motors 
and a-actinin cross-linkers, has shown that contractil- 
ity occurs above a threshold motor concentration (one 
myosin filament for every 30 actin filaments) and within 
a window of cross-link concentrations (90 — 270 a-actinin 
dimers per actin filament)^. It is somewhat unexpected 
that at high cross-link densities the bundled networks do 
not contract on the experimental timescale of an hour. 
This seems to imply a dramatic slow-down of cytoskele- 
tal kinetics due to the steric constraints associated with 
increasing crowding. 

The cytoskeleton in different organisms or even during 
different stages of the cell cycle exhibits a rich variety 
of viscoelastic propertie d 16 ' 17 . To accommodate both the 
adaptive behavior of dynamic remodeling, and the abil- 
ity to stabilize and be resistant to deformation, the cy- 
toskeleton should exhibit at least two mechanical phases: 
a plastic/fluid phase for internal reorganization of the 
cell, and an elastic/solid phase for mechanical support 
and tension transmission. 

In contrast to the strain-stiffening behavior— in re- 
sponse to a sustained stretching of cells or reconstituted 
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cross-linked actin gels, measurements on the responses 
to a transient stretching in the living cells^ demonstrate 
that the cell will fluidizc in response to the stretch, but 
will resolidify subsequent to this fluidization. This abil- 
ity seems to be insensitive to molecular details, and in- 
stead only depends on the proximity of the thermody- 
namic conditions of the cell to a solid-like state before the 
stretch. This observation implied that in addition to spe- 
cific signalling pathways, some response mechanisms are 
most likely controlled by non-specific actions of a slowly 
evolving network of physical forces. On the other hand, 
the universal behavior observed in the osmotically com- 
pressed cell^S highlighted the crowding-induced stiffening 
of the cytoplasm as the solid volume fraction is suffi- 
ciently high, and suggested an analogy to the colloidal 
glass transition. 

In an attempt to understand the interplay between the 
individual filament properties and network architecture, 
we study in this paper the equilibrium properties of a 
model cytoskeleton as an amorphous network of rigidly 
cross-linked nonlinear strings which also contains nodes 
with excluded volume. We incorporate the nonlinearity 
of interaction due to bucklable filaments into the force 
law, and characterize the cytoskeletal architecture by 
cross-link density and network connectivity. A mean- 
field level investigation, within the framework of the self- 
consistent phonon method and a density functional for- 
mulation, reveals a diversity of mechanical phases and 
a number of possible transitions in between which can 
be controlled by biophysical parameters. This phase dia- 
gram may shed light on our understanding of cytoskeletal 
remodeling in response to mechanical or chemical stim- 
uli. We also show the possibility of a glass transition 
in this model system and how the network connectivity 
modulates the transition densities. 



II. MODEL AND METHOD 
A. Model 

We model the F-actin bundles as individual nonlinear 
elastic strings; they are capable of resisting tensile forces 
by stretching (beyond the relaxed length L e ), but are un- 
stable under compression so that these forces may cause 
the strings to buckle. The actin-binding proteins bundle 
and crosslink F-actin to form an amorphous network. To 
account for the excluded volume of F-actin and the ABP 
aggregates, to the lowest order, we model the network as 
being crosslinked by hard "beads"; each bead serves as 
a compact rigid subunit that concentrates the volume of 
the F-actin and ABP aggregates centered on that bead. 
The network elasticity then comes from the remaining 
inter-bead F-actin molecules which are now taken to be 
volumcless. To higher order, the elasticity within the 
beads can also be included by introducing softness in 
the repulsion. In this way we effectively decompose the 
inter-bead interaction into purely excluded volume and 



nonlinear clastic contributions, and the model potential 
energy V(r) between a nearest-neighbor pair of beads is 
given by 

pV(r) =AQ(d-r) + ^(r- L e ) 2 0(r - L e ) , (1) 

where Q{x) is the Heaviside step function. The limit 
A — > oo indicates the hard-sphere (HS) repulsion, while 
7 measures the rigidity of the inter-bead F-actin. Tem- 
perature dependence enters this model only via the com- 
bination /?7 where (3 = I/UbT. Here d denotes the HS 
diameter of the beads and L e marks the onset of elastic- 
ity. Thus L e > d define a buckling regime (d < r < L e ) 
where no load is imposed. We show a sketch of the model 
system (Fig. [T^.) and a schematic of the nonlinear inter- 
action (Fig. [lb). 




FIG. 1: Illustrations of the model system and the nonlinear 
interaction, (a) The beads (blue spheres) interconnect the 
F-actin (red straight or green squiggly lines) into an amor- 
phous network. Red straight lines stand for tense/stretched 
bonds and green squiggly lines for loose/buckled bonds. The 
arbitrarily chosen central particle (yellow-stared) is connected 
with its nearest neighbors within the first shell of the radial 
distribution function (the dotted purple circle), (b) Nearest- 
neighbor interaction versus radial separation. 

We incorporate the filament nonlinearity to the extent 
that it effectively encodes the asymmetric response of 
the filaments to stretch and to compression; in this sense, 
our model is more "coarse-grained" than the well-studied 
"scmi-ficxible fiber" model^I; the latter treats the bend- 
ing degrees of freedom of F-actin based on the contour of 
the filaments using a continuum description. This bend- 
ing effect, stated in our language of inter-bead interac- 
tion, essentially introduces a finite resistance/restoring 
force to compression in our originally buckling regime. 
It can be easily incorporated, yet we do not expect any 
qualitative modification on the mean-field level predic- 
tions made by our current model. 

Similar buckling bonds were used to study the sta- 
tistical mechanics of a "cat's cradle" built on a regular 
lattice by two of us several years ago22. Here we also 
include a HS repulsion to account for the excluded vol- 
ume of the F-actin and ABP aggregates, and also assume 
the cytoskeleton adopts an aperiodic amorphous struc- 
ture characteristic of the disorder of biological cytoskele- 
tal networks. Thus in our model system, localization is 



3 



achieved not only at a high concentration of beads via 
hard-core repulsions (topological caging), but also will 
occur upon network expansion due to bond stretching. 
The interplay between the lioiiliiiearity of interaction due- 
to bucklable bonds and the network architecture renders 
the thermodynamic state diagram nontrivial in terms of 
diverse mechanical phases. These equilibrium phases in- 
clude persistent uniform liquid-like states, regions with 
a coexistence of frozen and liquid phases, and the possi- 
bility of a martensitic-like phase transition that signals a 
spontaneous symmetry breaking. 



The average length of F-actin in an in vitro network 
can be adjusted by controlling the concentration of cap- 
ping proteins (for example the concentration of gelsolin); 
a higher molar ratio of capping protein to actin results in 
shorter F-actin on average. Increasing the crosslink con- 
centration promotes the formation of F-actin and ABP 
aggregates^ 3 -, and in turn increases the bead density. The 
bead density determines the average spacing between the 
beads, or cquivalently, the average end-to-end distance of 
F-actin (r). When r < L e the elasticity of the fiber is 
cntropic in origin, whereas the intrinsic elastic modulus 
of the fiber dominates when r > L e . The effective stiff- 
ness of the interbead F-actin (/?7) is also variable; it is 
enhanced when filaments are bundled together to form 
structures with larger diameters^. Thermal fluctuations 
play a smaller role as filaments become stiffer. 



To get a feeling about how the relevant biophysical 
parameters work out in our model system, we can first 
estimate the bead density p from experimental actin con- 
centration. The typical actin concentration used for in 
vitro networks is 23.8/iM which, if wc assume the actin 
monomers to be spheres of diameter 5nm, corresponds to 
a volume fraction of 0.9 x 10~ 3 . The ABPs take up neg- 
ligibly small volume compared to that taken up by actin. 
We further take the bead size to be 1/xm and the F-actin 
as slender rods of cross section 5nm x 5nm, then if each 
bead concentrates five F-actin of length 6/im, p = 1.2 is 
needed to reach the given volume fraction; if five F-actin 
of length 10/mi per bead then p = 0.9 (taking the bead 
diameter d to be the length unit). Thus p is adjustable, 
ranging from 0.1 to 1.4, by varying actin concentration 
and/or crosslinking properties. Since the filament ag- 
gregates are not perfectly dense-packed within the bead, 
they only take up a portion of the assumed bead volume; 
the modeled hard-sphere repulsion between the beads 
may overestimate the excluded volume effect. Yet the 
qualitative phase behavior of this model system should 
not be altered. In a first approximation, the cytoskeleton 
determines the mechanical properties of a cell; since the 
elastic modulus of a cell is in the range of 10 3 Pa, the 
corresponding effective stiffness (3j is then estimated to 
be between 1 and 10 as converted into our model param- 
eters. (The characteristic bead size is d = lpm.) 



B. Mean field approximations and the self 
consistent phonon (SCP) method 

Instead of treating the bending degrees of freedom of 
F-actin directly, we focus on the motion of individual 
beads, located on the vertices of the network, in order to 
find the nonlinear elasticity as the bead density is varied. 
In the mean-field spirit, we tag a given bead as the cen- 
tral particle, and study its stability/response to the local 
mechanical environment. 

We build the model system on top of an amorphous 
structure. In contrast to a regular lattice that usually 
has a unique equilibrium configuration as well as a defi- 
nite coordination number for a specific lattice structure, 
amorphous systems must be appropriately averaged over 
non- vibrational disorder in the lattice which we take as 
quenched. A further mean-field approximation will be 
made in the present analysis to avoid detailing the con- 
figurational complexity of a random network. We will 
summarize the underlying topology of the amorphous 
solid in an assumed isotropic pair distribution function 
g(r) for the fiducial configurations of the system. Such 
a treatment has also been used for molecular structural 
glasses. We then define nearest neighbors as those beads 
that sit within the "first shell", i.e., up to the first mini- 
mum of the equilibrium radial distribution function g{r), 
and assume that interaction only exists between nearest- 
neighbor pairs. 




t 2 3 4 5 o 02 0.4 0.6 0.S I 12 

K P 

FIG. 2: (a) Radial distribution function g(R) vs the radial 
separation R for a series of particle density p. As density 
increases, oscillatory amplitude gets larger along with a phase 
shift toward the core, (b) Coordination number vs particle 
density. 



For liquids above the melting point, gns is well de- 
scribed by the Percus-Yevick approximation; the Verlet- 
Weis correction improves the behavior of g(r) near the 
core and dampens its oscillations at large r, giving accu- 
rate modifications especially at high densities^. We fol- 
low the procedures of Verlet and Weis, and note that co- 
ordination number in an amorphous structure depends on 
particle concentration, because the amplitude and phase 
of the oscillations in radial distribution vary with number 
density /packing fraction of the particles (See Fig. (2ja)). 
As shown in Fig. [21b), an increasing number of nearest 
neighbors are accommodated within the first shell when 
particles arc packed denser. This increase is almost linear 
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when p is low, and slows down as p approaches 0.6; the 
modest increase above densities of 0.8 is bounded by the 
random close packing value of the coordination number 
which is around 14. 

In the context of our model network, the amorphous 
topology effectively contains the physical aspects of bond 
breaking and/or ABP detachment upon network expan- 
sion. These aspects were absent in the lattice setting 
which was studied earlier: as the density of beads de- 
creases (or equivalently as the network expands), coor- 
dination number drops off resulting in a smaller number 
of interacting neighbors or an effectively weaker network 
connectivity, even though the bond connections between 
each nearest-neighbor pair are assumed to be permanent 
at each given bead density. 

As an approach of thermal stability analysis in an equi- 
librium system, the SCP theory was first developed to 
treat the anharmonic effects of hard-sphere crystals. The 
basic idea is to introduce a reference harmonic system, 
and then obtain the effective potential felt by each tagged 
particle, by averaging the pair interaction over the as- 
sumed Gaussian fluctuations from all its neighbors. This 
procedure should give back the assumed harmonic po- 
tential of the typical particle. The resulting coupled set 
of self-consistent equations allows an iterative scheme to 
determine the generally site-dependent force constants. 
Among several schemes, Fixman's SCP method based on 
a systematic expansion in Hermitc functions has proven 
an efficient and especially robust procedure^. Recently, 
this technique has been applied to network glasses^, and 
also by two of us to study motorized particle assembly^ 
to analyze the far-from-equilibrium dynamics like that of 
the living protoplasm using a local feedback scheme. 

We apply the SCP method to quantify the responses of 
various mechanical phases under varying physical condi- 
tions. At low concentration of beads, the effective attrac- 
tion due to stretched springs dominates, thus the tagged 
particles are localized by "bond trapping" ; whereas in the 
high concentration limit, HS repulsion dominates and 
results in a glassy /jammed state owing to "topological 
caging". To investigate both localized phases and the 
intermediate states that bridge the transition, we make 
a Gaussian local density profile ansatz with a single pa- 
rameter, i.e. we describe the time-averaged density con- 
figuration as a sum of Gaussians representing thermal 
vibrations of particles about the fiducial sites 
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In the present work, the force constants (inverse mean 
squared displacements or localization strengths) will 
all be taken to be equal, but it is not difficult to allow 
spatial variation^, i.e., dependence on the index i. 

In the independent-oscillator version of the SCP the- 
ory, the effective potential between two interacting par- 



ticles is given by 
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which may be explicitly written as 
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Here, R denotes the averaged equilibrium separation be- 
tween interacting particles, and a represents the homo- 
geneous localization strength. 

Taylor expansion of the effective interaction up to the 
second order gives a self-consistent relation for a: 

a=^[ d 3 Rg(p, R) Tr[W(3V e "(R, a)] , (5) 
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J dRR 2 g{ Pl R)V 2 l3V en (R,a) . (6) 

Here R* marks the position of the first minimum of the 
radial distribution function g(R). The HS diameter d, as 
the lower limit, is taken to be the unit of length. 



C. Thermodynamic ramifications of our model 
system 

Wc shall show in a moment that in our model system — 
a rigidly cross-linked nonlinear-elastic network — there 
exist at least two different localized phases: one mod- 
estly localized phase describes the weakly trapped mo- 
tion due to bond stretching. We refer to this state as 
the "liquid-like state", in view of its considerable mo- 
bility and small localization strength (au q ). The other 
more strongly localized state corresponds to the jammed 
motion within topological cages of neighbors. This solu- 
tion depicts a "glassy state" exhibiting highly restricted 
vibrations (a g i). 

Traditionally, microscopic treatments of liquids take 
the view that since the liquid structure is dominated by 
repulsive forces, it is desirable to develop perturbation 
theories based on a HS reference system and then find the 
optimal parameters for it. In this spirit, the Hclmholtz 
free energy can be obtained by adding a first-order per- 
turbation to the free energy of the corresponding HS sys- 
tem. On the other hand, the Carnahan-Starling equation 
of state gives accurate values of the reference free energy 
at moderately high densities. In sum, the free energy for 
the liquid-like state is given by 

f lig = ^ = (inpA 3 - 1) + f (Zestf) ~ 1) 
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Here the first term gives the entropic cost when all the 
nearest-neighbor pairs are bonded, with A denoting the 
thermal wavelength. The second term is the excess free 
energy of the HS reference system where the compress- 
ibility factor Z is given by Carnahan-Starling(CS) equa- 
tion of state. The last term involves the energetic contri- 
butions from the attraction due to bond stretching; note 
the HS part is carefully deducted. 

We use SCP theory to describe the free energy of glassy 
configurations^^!. The expression we use for an individ- 
ual glassy configuration is 



The number of stretched bonds can then be obtained by 
averaging over non-vibrational disorders, i.e. configura- 
tional degrees of freedom: 
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Here the first integral gives 
ergy after double averaging 
ations and over topological 
ond term (inside the curly 
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over thermal fluctu- 
disorder. The sec- 
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which accounts 
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for the effect of cell constraint; Z p represents the mo- 
mentum part of the partition function, and we choose a 
cubic cell with side length D = p~ 1 / 3 /2 for convenience. 
Finally, Sf (taken to be 0.224) is a numerical correction 
to the entropy of HS crystal system near face-centered- 
cubic (fee) close packing obtained via SCP approxima- 
tion, which then gives extremely accurate free energies 
near melting. 

The pressure can be evaluated by numerically differen- 
tiating the liquid free energy: 



P = P 



r, JUq 
op / T N 



knT. 
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Random first order transition theory identifies the con- 
figurational entropy with the difference between the free 
energy of the highly localized glass solution and the liq- 
uid, i.e. 
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Another interesting quantity that can be evaluated 
within the SCP theory is the number of force-bearing 
bonds, i.e., those that have a length exceeding the elas- 
ticity onset L e . In three dimensions, the structure- 
dependent probability for a single bond to be elasti- 
cally stretched (or equivalently, the fraction of stretched 
bonds, in the mean field context) is given by 
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q3(a,p)=P dRg(p,R)q 3 (a,R). (13) 

J 1st shell 

This double integral accounts both for the contributions 
from the fluctuation of individual bonds (small a indi- 
cates strong fluctuation) and from fluctuations in the 
underlying topological structure (short L e enlarges the 
radial range that contains stretched "fiducial bonds" ) . 



III. NUMERICAL INVESTIGATIONS 

A. Localization strength 

The fundamental quantity that characterizes the di- 
verse mechanical phases in our model is the localization 
strength, or the force constant of the emergent Einstein 
oscillator, a. We may plot a against bead density p 
and/or effective stiffness /?7 of the F-actin. We mea- 
sure lengths and energies in units of d and (3, then the 
corresponding dimensionlcss quantities are taken to be 
p* = pd 3 = p and 7* = fijd 2 = f3j. We start with sev- 
eral representative one dimensional plots that come from 
vertical slices of the two dimensional a-surface, and we 
first focus on the liquid-like solution which is absent in 
the pure HS system. 




1 1.2 1.4 



(10) FIG. 3: Liquidlike localization strength au q versus 
bead density p. (a) For 7* = 2 with L e = 
1.2 (blue), 1.5 (green), 1.8 (red); (b) for L e = 1.8 with 
7* = 2 (red), 10 (green), 30 (blue). 

Referring to Fig. [2ja), we observe that the first shell 
(R*) of the radial distribution shrinks from 1.95 to 1.3 
as the bead density increases from 0.2 to 1.2. When 
L e < R* , the separation (R* — L e ) determines how many 
nearest-neighbor fiducial sites are found beyond L e and 
on-average have tense bonds; in this case, both the un- 
derlying topology and thermal fluctuations contribute to 
the localization. On the other hand, however, if L e > R* , 
all fiducial sites of nearest neighbors fall inside the sphere 
of radius L e and on-average result in buckled strings; in 
this case fluctuations are the only source of stretching and 
thus of localization. This information is also encoded in 
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the threshold density p t h beyond which a stable au q so- 
lution no longer exists. In sum, longer onset length leads 
to weaker overall localization as well as lower threshold 
density, as shown in Fig. Ela). A similar effect is pro- 
duced by low rigidity 7 as can be seen in Fig. (3jb); since 
smaller 7 indicates a broader and shallower confining well 
in which particles are more loosely tethered thus being 
less localized, and the corresponding liquid-like solution 
becomes unstable at lower p t h- 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 
P 

FIG. 4: Liquidlike localization strength versus bead density 
for L e = 1.2 with /J7 = 2, 10,30. Liquid-like solution is still 
distinct at ftj — 2, whereas a rapid crossover to glassy behav- 
ior occurs for /J7 = 10, 30. 

The compromise between the number of contributing 
neighbors and the degree of stretching produces a non- 
monotonic density dependence of au q at a given L e , and 
the most efficient localization is achieved at an interme- 
diate density around 0.6. At p < 0.6, the individual 
bonds become less likely to be stretched as the beads 
pack more densely, whereas the total number of bonded 
neighbors increases much faster with increasing density, 
the combined effect thus leads to a stronger localization 
as density grows. On the other side when p > 0.6, the 
number of bonded neighbors almost saturates, while the 
thinner first shell of g(r) indicates less tightly stretch- 
ing or deeper buckling as density increases, whereby the 
localization strength decays accordingly. 

Note that this non-monotonic density dependence only 
occurs for a sufficiently long L e or in the low-,07 regime 
for a short L e . For sufficiently high ^7 liquid- like solu- 
tions are no longer distinct, instead a rapid crossover to 
high-a solutions is observed, as is exhibited in Fig. 0] 



with bead density for various L e as long as stable liquid- 
like solution exists (Fig. EJa)); at high 7*, however, close 
to the most localized regime for individual beads, even 
though the coordination number still gently increases, 
the weak fluctuation (high a) strongly suppresses bond 
stretching, thus leading to an intermediate decline of q% in 
its p-dependence, as seen in 7* = 10, 30 cases (Fig.EKb)). 




p p 



FIG. 5: Average number of stretched bonds 53 vs bead density 
p. (a) 7* = 2 for L e = 1.2,1.5,1.8 (top to bottom); (b) 
L e = 1.8 for 7* = 2, 10, 30 (top to bottom). 



2. free energy profile — consistency with SCP theory 



free energy profile (L e =1.8, Py=30) 
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FIG. 6: Free energy fu q versus the order parameter — 
localization strength a with L e = 1.8, ^7 = 30 for p = 0.2-1.2 
(bottom to top). As the bead density p increases, the overall 
profile shifts upward with the valley indicating the equilib- 
rium solution auq. 



B. Thermodynamics 

1. q-i — role of localization 

Recall that explicitly counts the average number of 
tensely bonded neighbors. Its dependence on elasticity 
comes from the self-consistcntly determined localization 
strength a and the intrinsic cutoff L e for a bond to be 
stretched. For a given p, more elastically-bonded neigh- 
bors become available as L e drops (Fig. [3(a), bottom to 
top) and/or 7* decreases (Fig. 0(b), bottom to top). For 
a moderate 7*, the number of stretched bonds increases 



As can be seen when we compare Fig. [5] with the cor- 
responding curves in Fig. [3jb), the minima on the free 
energy profile F(a) for the liquid-like state coincide with 
the corresponding auq solutions obtained via the SCP 
method. As the bead density increases, the overall profile 
shifts upward. At a given density (p = 0.6 here), rais- 
ing L e (Fig. Efa)) or lowering 7* (Fig. [3(b)) causes the 
liquid-like state to be increasingly favorable (deeper and 
sharper valley) with lower localization strength (bottom 
of valley being shifted leftward). As expected, 7* only 
markedly affects the low-a regime (a < 10), since when 
particles are sufficiently localized (due to short L e or high 
p), bond stiffness only plays a minor role in altering the 
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(a) 



free energy profile (fr(=2. p=0.6j 



(b) 




free energy profile tI. r =!X. p=fW) 




FIG. 7: Free energy versus localization strength, (a) 7* = 2 
for L e = 1.2, 1.5, 1.8 (top to bottom); (b) L e = 1.8 for 7* = 
2, 10,30 (bottom to top). 




interaction strength. In contrast, varying L e not only 
modulates the low-a regime, but causes a nearly uniform 
upshift of the high-a portion of F(a); this observation is 
consistent with our previous statement that when ther- 
mal fluctuations are weak (high a) topology dominates; 
shorter elasticity onset L e indicates more stretched bonds 
at a given density and thus enhanced interaction. 



3. pressure 

We plot pressure versus bead density with 7* = 2 for a 
series of onset lengths in Fig. [8j At the low-density end, 
most of the bonds are stretched and the system tends to 
shrink, resulting in a very low pressure. At the oppo- 
site end, in a densely packed system HS repulsion takes 
over the major role; the steeply increasing pressure re- 
flects the rising difficulty in rearrangement and the con- 
sequent soaring resistance to compression. In contrast 
with the nearly universal behavior at both extremal-ends, 
parameter-sensitive features emerge in the intermediate 
density regime, where shorter L e enhances the effective 
attraction and thus lowers the pressure of the system. In 
particular, a non-monotonic behavior is observed to oc- 
cur for pliable springs (low 7*) with early onset (short 
L e ), as is exhibited for the case of L e = 1.2. 

It is natural to expect the possibility of a negative pres- 
sure at a modest bead concentration in view of the ten- 
dency of the network to shrink due to attraction. We plot 
in the inset of Fig.|8jthe contribution of attraction (p a ttr) 
to the total pressure (ptot) for several onset lengths, and 
Pattr docs indeed exhibit a negative valley at intermedi- 
ate concentrations. It is clearly seen that a short onset 
length is necessary for a notable contribution from attrac- 
tion and thus for the emergence of this non-monotonic 
behavior of the total pressure. 

On the other hand, the negative contribution owing to 
attraction is eventually overwhelmed by the positive part 
due to a more rapid increase in repulsion as p increases, 
and the total pressure (ptot) remains positive over the 
whole density range. This behavior is originated from 
the specific form of the free energy functional that we 
have used for the model system, where the correction due 



FIG. 8: Pressure versus bead concentration for soft strings 
(/?7 = 2) with various onset lengths L e = 1.2, 1.3, 1.4, 1.5 
(bottom to top). Main plot: total pressure; inset: pressure 
due to attraction. 



to attraction enters as a perturbation to the dominating 
repulsive part. This result might require modification to 
better accommodate the low-concentration regime where 
attraction becomes important. Another reason for this 
behavior may be related to the g(r) we have used for 
the radial distribution which drops sharply when move 
away from the core; bonding effects might lead to a fat- 
ter tail of g(r) and an enhanced contribution from at- 
traction. Nevertheless, we feel these insufficiencies of the 
approximations should not modify the mechanical and 
thermodynamic properties qualitatively. 

The thermodynamic considerations discussed above 
help with understanding of the physics underlying the 
predicted phase behavior in terms of the order parame- 
ter a. We will discuss these in detail next. 



C. State diagram 

1. Overview 

In our model network, we assume permanent bond con- 
nections between the nearest-neighbor pairs. The inter- 
play between these intrinsic constraints and the thermal 
fluctuations gives rise to an inhomogencous distribution 
of tensile forces throughout the network. Such force het- 
erogeneity exhibits as a (uniform) non-vanishing tether- 
ing strength in our mean field context. This explains the 
prediction that a completely ergodic fluid phase, which 
is allowed for the pure hard-sphere system with short- 
range attraction and that has strictly diffusive behavior 
for long times i.e. a — > as t — > 00, never occurs for the 
equilibrium network structure; the lack of utter freedom 
in locomotion leads to a finite a over the whole span of 
bead concentrations. 

There exist two characteristic densities in our model 
system: the threshold density p t h above which the ho- 
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mogcncous liquid-like solution is no longer stable, and 
the critical density p cr which signals the emergence of a 
glassy state. These densities help define the boundaries 
between diverse mechanical phases: 

• P < Pth, Per'- only an q exists, which describes the 
liquid-like loosely tethered phase; 

• Pth < P < Per- the mean-field otu q is no longer sta- 
ble, and the glassy state has not yet occurred; bifurcation 
to a low-a solution may occur, exhibiting a symmetry 
broken phase which we shall call the "martensitic-like" 
(ML) phase; 

• Per < P < Pth- both ctuq and a g i exist, depicting 
the transition state with presumably coexisting phases 
that implies a macroscopic number of configurational de- 
grees of freedom for structural rearrangements. 

• p > pth, Per'- the repulsive-glass phase dominates. 



2. Features of the various phases 

Our SCP calculation has found five distinct phases 
in our model system: the liquid-like (LL) phase, the 
crossover (CO) phase, the repulsive-glass (RG) phase, 
the multiple-solution (MS) phase, and a martcnsitic-likc 
(ML) phase. 

the liquid-like (LL) phase and crossover (CO) phase 

We choose two particular onset lengths based on their 
position with respect to the boundary of the nearest- 
neighbor shell. 

For L e = 1.8, this cutoff stands outside the nearest- 
neighbor shell for p : 0.2-1.2, which means that almost 
all the "fiducial bonds", i.e. bonds that connect fidu- 
cial sites of nearest neighbors, are buckled, and become 
more floppy as bead density grows; in this situation once 
the coordination number saturates at p ~ 0.6, increas- 
ingly deeper buckling results in a decrease in au q as p 
increases. Owing to the overall buckling of the bonds, 
higher bond stiffness and thus smaller fluctuations lead 
to a smaller probability for the bonds to tense up, and 
an even steeper drop in au q occurs as density increases 
(Fig. EJb)). Therefore, the distinct liquid- like solutions 
persist over the whole 7*-range of interest (as long as 
P < Pth) for the case of a long onset length. 

For the case of L e = 1.2, however, since this cut- 
off remains inside the nearest-neighbor shell all the way 
through p : 0.2-1.2, there always exists a fraction of fidu- 
cial bonds being stretched beyond their relaxed length. 
In this situation the descending branch of au q at inter- 
mediate densities only occurs for very soft springs, i.e. in 
the I0W-7* regime. Furthermore, taking advantage of the 
persistent fraction of stretched (fiducial) bonds, sufficient 
stiffness beyond a threshold value (marked by the phase 
boundary p co ) would help enhance localization, and fa- 
cilitates a smooth crossover from elasticity-dependent 
liquid-like behavior to geometry-dominant glassy behav- 
ior (as seen in Fig. HI) ; we shall refer to the states showing 
such behavior as being in a "crossover" (CO) phase. 

the repulsive-glass (RG) phase 



Current simulations performed on increasingly longer 
timescales make it possible to compare the equilibrium 
properties based on the present theoretical predictions 
with simulation outcomes in the long-time limit. To in- 
vestigate the long-term fate of attractive glasses, simu- 
lations of glassy arrest in hard-core particles with short- 
range attraction 3 ^ were performed over a waiting-time- 
independent window of up to 10 6 MD units, which is a 
few orders of magnitude longer than those reached by 
previous experiments or simulations (~ 10 3 MD units). 
They found that even if the short-range attraction gener- 
ates a transient plateau in the time-evolution of the cal- 
culated mean-squared displacement (i.e. inverse a) owing 
to breaking and reforming of nearest-neighbor "bonds" , 
the long-time behavior of "bonded" and "nonbonded" re- 
pulsive glasses converges, suggesting that in the long run 
particles are ultimately confined by their topological cage 
of neighbors. 

Our equilibrium calculation on a long-rangc-attractive 
HS system consistently shows that as long as the (nearly- 
universal) critical density is reached, there would emerge 
the repulsive-glass behavior being almost independent of 
the attractive strength 7*. Such independence is found 
for various onset lengths. 

the multiple-solution (MS) phase 

The simultaneous presence of distinct au q and a g i so- 
lutions, according to the SCP analysis, implies a non- 
vanishing configurational entropy in this phase region. 
This further signifies the availability of configurational 
degrees of freedom for structural rearrangement. As 
shown in previous examples, a longer onset length (L e = 
1.8) allows for distinct au q solutions over a larger (3j 
range and thus allows the system to explore more en- 
ergetically favorable configurations; on the other hand, 
the emergence of crossover behavior at moderate /?7 and 
the consequent shrinkage of the MS regime induced by 
early onset of elasticity (L e = 1.2) suggests a rapid loss 
of configurational entropy as bonds become stiffcr. 

the martensitic-like (ML) phase 

As verified by our SCP calculation, a completely er- 
godic fluid phase, presented by a sticky HS system with 
weak attraction at a low density, does not occur for an 
clastically bonded HS system, so that a is always fi- 
nite. Yet finite localization strength does not ensure a 
homogeneous structure. In our model system, with a 
large elasticity onset L e (compared to the lattice spac- 
ing) and a low bond stiffness 7* , bifurcation in au q takes 
place before the emergence of glassy behavior (i.e. below 
p cr ) and may allow the existence of a spatial-symmetry- 
broken phase characterized by a stable pair of liquid- 
like solutions. The occurrence of bifurcation was also 
found in an earlier study by Shen and Wolynes with a 
pure "cat's cradle" built on a regular lattice. The pos- 
sibility of self-generated spatial heterogeneity associated 
with such a mechanical instability makes it interesting to 
study what kind of structural phase such a destabilized 
system would actually develop into, and to quantify the 
associated conformational changes. The resultant me- 
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chanical/structural phase might be related to the orien- 
tational order often observed in F-actin networks of the 
cytoskclcton. 

We performed explicit molecular dynamics (MD) 
simulations on a body-centered-cubic (BCC) net- 
work of n 3 unit cells with pure nonlinear elas- 
ticity, as described by the Hamiltonian H = 
5 E< E< 3 > b (\n - r!\ - Lef e(\Tt-¥l\-L e ). 
Here i,j label the nodes and (• • •) represents sum over 
nearest neighbors; 7 denotes bond stiffness and L e 
elasticity onset as before. The bipartite nature of the 
BCC structure allows us to divide the original BCC 
network into two interpenetrating simple-cubic (SC) 
subnetworks (see an illustration of the subnetwork divi- 
sion in Fig. IHta)) , so that each node on one subnetwork 
interacts with its 8 neighbors on the other subnetwork, 
and the neighbor list never changes. 

The quality of the motion — whether it is oscillatory 
or monotonic — depends on the relative contribution of 
the inertial forces (that tend to produce oscillations) 
and the viscous forces (that tend to damp the oscilla- 
tions out). It turns out that inertial forces are usu- 
ally very small at the microscopic and molecular lev- 
els, so that the overdamped limit usually applies^. We 
thus carried out the simulations with stochastic dynam- 
ics in the overdamped limit as described by the Langevin 
equation dfi/dt = (l/T)/j(i) + rfi(t). Here / represents 
the deterministic force due to the nonlinear elastic in- 
teraction, and the force exerted by the fluid particles 
divides into two parts: the average viscous force — Til 
and a random force £(t) = Tfftt) whose time average is 
zero. We assume a Gaussian white noise that satisfies 
= 2D5ij5(t - t'). As usual, V and D denote 
the drag coefficient and diffusion constant, respectively. 
Lengths are expressed in units of L e . We applied periodic 
boundary conditions. 

We first show the evolution of the mean squared dis- 
placement (MSD) with respect to the initial equilibrium 
positions of the nodes. We simulated a system of size 
5 3 with soft bonds (b = /?7 = 1) at a low temperature 
(f3 = 30). When the elasticity onset L e is compara- 
ble with the initial mesh size, i?, of both subnetworks, 
clastic stretching is immediately felt as soon as thermal 
buffeting displaces any of the nodes from their equilib- 
rium locations, and therefore, the whole network just 
wiggles about the underlying BCC lattice, which defines 
the unique minimum of the system's energy landscape. 
As shown in Fig. M)>) , hi this case the two subnetworks 
act in a concerted manner with their mean squared dis- 
placements fluctuating almost "in phase" , suggesting the 
whole BCC network does not show symmetry breaking. 

When the system parameters are modified to be in 
the regime of bifurcation by increasing the L e /R ratio 
to 10 (other parameters unchanged), however, the MSD 
exhibits a much larger average amplitude (after a steady 
value is reached) and stronger fluctuations around it; fur- 
thermore, the two subnetworks fluctuate in a correlated 
manner but are almost completely out of phase with each 
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FIG. 9: (a) An illustration of the subnetwork division. The 
bipartite nature of the original BCC lattice allows for sep- 
aration of one SC subnetwork (red) from the other (blue). 
Each node on one subnetwork (blue sphere) interacts with its 
8 neighboring nodes (red spheres) on the other subnetwork, 
(b) The MSD of both subnetworks (n — 5) versus simulation 
time in MD units for a moderate elasticity onset to mean 
separation ratio L e /R = 1. T = 50, P = 30. 
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FIG. 10: Effect of the system size on the behavior of MSD in 
the regime of bifurcation. Both: L e /R — 10, F — 50, P = 30. 
(a) n = 5; (b) n = 10. 



other, see Fig. [TUf aV In this case, the initial network 
presents no elastic constraints on the nodes, thus the sys- 
tem is free to expand until it reaches a steady size where 
the mean separation between the bonded neighbors is 
comparable with the length for elasticity onset; the cor- 
related fluctuations result from alternate distortions be- 
tween the subnetworks within a system of moderate size 
as 5 3 . As we enlarge the system to n = 10 (Fig. fTUTb)), 
the steady amplitude of MSD maintains (V M SD ~ 0.6 
stretching could occur under thermal driving), yet the 
fluctuations about the average become much weaker. 

The above observation raises the intriguing possibility 
of a "martensitic-like" phase consisting of (frustrated) 
domains with (complementary) distortions. To exam- 
ine further this possibility, we also queried vectorial in- 
formation about the directions of motion (lost in MSD) 
by examining the evolution of displacement vectors. We 
divided the simulation cell evenly along each initial di- 
mension to get 8 domains each with the same number 
of nodes, and traced the subsequent motion within each 
domain. We chose a time window t^D = 4000-6000 
after the steady size was reached, and used the run- 
ning time averaged position of the nodes as the reference 
with respect to which the displacements were defined. 
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We present in Fig. [TT] the bulk-averaged (left column) 
and domain-averaged (right column) displacement com- 
ponents. The evolution of the displacement vectors in 
these two cases shows several contrasting features. First, 
the maximal amplitude of the domain-averaged displace- 
ments is about one order of magnitude larger than the 
maximal amplitude of the bulk-averaged displacements. 
Second, the bulk-averaged displacements of the two sub- 
networks are exactly oriented in the opposition directions 
with essentially the same amplitude, whereas in individ- 
ual domains two subnetworks move almost in the same 
direction yet with different amplitudes. These two con- 
trasting features support the picture of localized distor- 
tions with different orientations in different domains. In 
addition, the peak value of the displacement amplitude 
occurs at different times for different domains, which sug- 
gests "frustration" between the domains. 
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FIG. 11: The bulk-averaged and domain-averaged displace- 
ment components for both subnetworks, n = 10, L e /R = 
10, r = 50, /3 = 30. (a)-(c) average over the whole simulated 
region; (d)-(f) average over one of the 8 domains. (a)(d)x 
component; (b)(e) y component; (c)(f) z component. 



To roughly estimate the domain size and to quan- 
tify the orientational correlation within the domains, we 
made equal-width shells centered at each node (the di- 
ameter of the outermost shell was taken to be equal 
to the box size of the simulated system), and com- 
puted two measures: (1) the average projection of di- 
rector (i.e. the unit vector of the corresponding displace- 
ment) in each shell onto that of the central node, and 
then averaged over all possible central nodes. In mathe- 
matical terms, we monitor P m (t) = (1/N) ' 



(J-T 1 



(t) ) , where N is the total number of 



nodes and N m the number of nodes within the m th shell 
of the central node; f i and fj m are instantaneous directors 
of the central and in-shell nodes, respectively; (2) the av- 
erage tensor product of directors constructed as Q m (t) = 

(i/N) e£i m) ■ y) (it £l m =i (y ■ ^ (*))) • where v 

is the unit vector pointing from the instantaneous po- 
sition of node i to that of node j m . This measure thus 
reflects the degree of alignment of two directors along the 
line connecting their positions. 
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FIG. 12: The evolution of two measures of the orientational 
correlation for both subnetworks, n = 10, L e /R = 10, T = 
50, /3 = 30. (upper) the average director projection P for 2 
shells (a) and 4 shells (b); (lower) the average tensor product 
of directors Q for 2 shells (c) and 4 shells (d). In each panel, 
shown from top to bottom are values of the measure from the 
inner to the outer shell(s); red/magenta denotes subnetwork 
1, and blue/cyan denotes subnetwork 2. 

We show in Fig. [TJ] the evolution of these two mea- 
sures for the cases of 2 shells (left column) and 4 shells 
(right column). The average director projection P, in 
each shell, fluctuates about a steady value, which decays 
from around 0.2 to 0.05 as we go from the inner to the 
outer shell in the case of 2 shells (panel (a)). This spa- 
tial decay in average P indicates that the domain size is 
roughly half of the box size, i.e., 5 times the initial lat- 
tice spacing. As expected, the fluctuations in P become 
weaker for outer shells due to the larger number of nodes 
to be averaged over; this trend is clearly exhibited in the 
case of 4 shells (panel (b)). 

In contrast to the positive steady value of P in all 
shells, the average tensor product of directors Q becomes 
negative for the outer shell(s) (lower panels of Fig. IT2"j) . 
Since Q encodes the angular location of the in-shell node 
relative to the central node, a negative value of Q implies 
the two directors point in the opposite directions along 
the line that connects them. A stronger opposite align- 
ment would give a more negative Q value. Q thus serves 
as an indicator of the loss in orientational correlation. 
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Also, negative Qs arc smaller in amplitude than positive 
Qs. 

To make the situation more apparent, we visualized 
the ensemble-averaged displacement vectors at different 
angles and different distances from the central node for 
the case of 2 shells (Fig. [T3|) . We took the ensemble- 
averaged direction of motion of the central node as refer- 
ence. Within the first shell (colored in pink) the average 
displacement in any angular range gives a positive pro- 
jection on the reference direction, in other words, the 
displacement field is anisotropic and the movements oc- 
cur mainly along the same direction as that of the central 
node. The aligned movement is most significant near the 
"equator" . The displacement vectors in the second shell 
(colored in blue), however, no longer exhibit a preferred 
orientation, instead they reorient considerably. Close to 
the "south pole" , the outward movements are almost per- 
pendicular to the reference direction, whereas small in- 
ward movements occur near the "north pole" . The spa- 
tial decay in oricntational correlation is thus apparently 
observable and consistent with the features exhibited by 
our correlation measures. 




FIG. 13: Average displacement vectors (purple arrows) at dif- 
ferent angles and different distances (inner and outer shells) 
from the central node (red sphere). Shown for one cross sec- 
tion of the spherical region. 

In sum, in addition to the multiple solution regime 
which may be related to the saddle point solution 
pictured as ergodic droplets formed against a glassy 
background in a finite-range system, our mean-field 
level calculation also predicts a parameter regime where 
nonlincar-clasticity-induccd spatial inhomogeneity is ex- 
hibited through a "martcnsitic-like" phase with local ori- 
ented distortions. 



3. Typical examples 

We plot the 2D surfaces of au g and a g i against p*(=p) 
and 7* (=Pj) at a given L e to examine how these physical 
parameters modulate the phase boundaries. The contour 
maps are also projected on the bottom as reference for 
the upcoming state diagrams. We choose two particular 
values of the elasticity onset L e that characterize typical 
tense networks (L e = 1.2) and floppy networks (L e — 
1.8). 




FIG. 14: 2D surface of localization strength over the param- 
eter space extended by bead density (p*) and elastic stiffness 
(7*) in the case of L e — 1.8 characterized by mostly floppy 
bonds. Here p* runs through 0.1 to 1.4 and 7* ranges from 1 
to 30. We show the a surface and its contour map for liquid- 
like (a) and glassy (b) solutions. The color scheme indicates 
the relative measure of the a values; the highest value within 
a given range is colored as bright yellow and the lowest as 
black. 



We start with L e = 1.8 case. For this relatively long 
onset length, liquid-like and glassy solutions are quite 
distinct over the whole parameter space. As can be 
seen from Fig.QjJa), the liquid-like localization strength 
takes on a "hump" shape along the p-axis peaking around 
p = 0.6, and elevates smoothly in the 7* direction. How- 
ever, the stable au q solution terminates at a sharp bound- 
ary defined by pth{"f* , L e ) beyond which the mean field 
auq solution becomes destabilized. This instability re- 
gion, located in the high-/? low-7* corner, presents a 
quasi-triangular shape which indicates the increased need 
of stiffer bonds to stabilize the loosely arrested state as 
density increases. The ML phase may arise as a possi- 
ble consequence of this (clastic-)nonlincarity-induced me- 
chanical instability. On the other hand, a g i emerges at 
p cr ~ 1. The values of the critical density and of a g i 
are nearly independent of 7*, as shown in Fig. I14f bh 
such stiffness-independence arises from takeover of the 
dominant role by HS repulsion in reconfiguring a densely 
packed system. 

The corresponding state diagram on p-7* plane is dis- 
played in Fig. I15f a). The trajectories marked by pth 
(blue curve) and p cr (dashed line) unambiguously divide 
the state space into four distinct phase regions: distinct 
liquid- like (p < p t h,Pcr) and glassy (p > p th , p cr ) phases 
locate at the opposite corners diagonally, while the rest 
of the space naturally divides into the ML and the MS 
phases depending on whether 7* < 7* (p t h < Per) or 
7* > 7c (pth > Per), respectively. For a given L e , p th in- 
creases with 7* whereas p cr is topologically determined, 
consequently when 7* < 7* and pth < P < Per, the re- 
gion of the ML phase narrows down as 7* increases due 
to decreasing {p cr — p t h) until it disappears at 7* = 7* 
(pth = Per); when 7* > 7* and p cr < p < p th , the 
MS phase takes over and broadens as 7* rises because of 
growing (p th - p cr ). 

In addition to the four types of phases exhibited by 
the floppy network (L e = 1.8), the state diagram for the 
tense network (L e — 1.2) presents a novel phase bound- 
ary separating out a large region featuring a crossover 
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ness increase). Moreover, the ML phase doesn't emerge 
for tense networks since the instability is avoided by early 
elasticity onset. 

Beyond the upper boundary of this stripe-shaped 
phase region, the a-surface gently mounts up toward the 
high-7* high-p* direction and smoothly crossovers from 
the elasticity sensitive behavior to the glass-like behav- 
ior as p cr is approached; consequently, in this crossover 
regime p cr no longer marks a clear transition boundary. 
The resultant state diagram is presented in Fig. nST bV 




FIG. 15: The state diagrams of a typical floppy network with 
L e — 1.8(a) and a typical tense network with L e = 1.2(b). 
These diagrams are constructed against the corresponding 
contour maps of liquid-like and glassy solutions, summariz- 
ing all possible phases partitioned by transition boundaries. 
In the L e — 1.8 diagram, 7* marks the crossing point of 
pth curve (blue) and p cr (or pa) line (black dashed). The 
"mechanical" diagrams are further integrated with thermo- 
dynamic characteristics — laboratory glass transition density 
pG (green dotted line) and Kauzmann density pk (red full 
line). The floppy network exhibits higher pa and pK than 
the tense network does. 



behavior. When wc compare the 2D surface for au q 
with that for a g i over the whole parameter plane, we 
find that they almost coincide except for a stripe-shaped 
region in the low-7* high-p* corner. When we zoom in on 
this region (contour map shown in Fig. [ToTd)) that peels 
off the smoothly ascending a-surface (see upper panels 
of Fig. [16]) . we observe that the diagram appears like a 
"squeezed version" of that for the L e = 1.8 case; namely, 
auq also exhibits a non-monotonic density dependence, 
and pth locates the stability limit of such mean-field so- 
lution at each 7*. As for the glassy state, however, a g i 
proceeds with its 7*-independent behavior as soon as the 
density exceeds p cr ~ 1 (see Fig.lWc)). Therefore, simi- 
lar to what was seen for the case of L e = 1.8, this stripe- 
shaped region presents LL, MS and RG phases. Note- 
worthy is that short L e dramatically reduces the region 
corresponding to the MS phase, indicating a rapid loss of 
the configurational degrees of freedom (upon bond stiff- 



FIG. 16: 2D surface and contour map of the localization 
strength for the case of L e — 1.2 that features a persistent 
fraction of tense bonds. Upper panels: liquid-like solutions 
up to auq — 100(a) and au q — 350(b). Both the peeling- 
off stripe-shaped region for distinct liquid-like solutions and 
the smooth crossover to glassy behavior are explicitly dis- 
played, (d) Amplified contour map for the stripe-shaped re- 
gion; (c) glassy solutions over the whole parameter regime. 



4- Characteristic densities 

Due to the bonding constraints inherent in a network 
structure, the stabilized state in our model always pos- 
sesses a finite localization strength, i.e. a > 0. Thus in 
this case the dynamical transition density p cr (or called 
Pa as in literature) , rather than being the lowest density 
to give a non-zero a as occurred in pure or sticky HS sys- 
tems, is defined as the lowest density to trigger persistent 
high-ct solutions over the whole 7* range of interest. In 
our model system, p cr signals the emergence of an ex- 
tensive number of glassy metastable states, yet does not 
mark the termination of bonding effect. 

While the SCP theory alone allows us to find pa, the 
"Kauzmann density", px 1 at which the Hclmholtz free 
energies of the liquid-like and glassy phases match and 
the configurational entropy ceases to be extensive also re- 
lies on the specific form of the free energy functionals we 
use for these two types of states. The ratio pa/ Pk dimcn- 
sionlcssly characterizes the thermodynamic aspects. To 
connect to the kinetic laboratory glass transition, we note 
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the laboratory transition is defined to occur when the vis- 
cosity reaches 10 4 Poise. Random first order transition 
(RFOT) theory predicts this to be when the configura- 
tional entropy is about 1.0 fcg per particle. To translate 
our thermodynamic results to the laboratory transition 
density, we will therefore mean by pa the density where 
the liquid and glass free energies differ by l.OfcsT per 
particle. Despite the universality of configurational en- 
tropy at laboratory transition experimentally confirmed 
in a wide variety of molecular glasses, this universality 
must be examined further to see if it is valid in the cy- 
toskeletal system which is an active biological material. 
Here we just use this fiducial entropy to discuss the qual- 
itative features of our model system. 

5. Possible transitions 

Visual inspection shows that the phase partition in the 
upper portion (7* > 7*) of L e = 1.8 state diagram ex- 
hibits qualitatively identical behavior to that found near 
the bottom of the stripe-shaped phase region in the di- 
agram of L e = 1.2 case, again indicating that high 7* 
and short L e are comparably competent for making ef- 
fectively more tense bonds. In this regime, melting from 
RG via MS region to LL state is expected as p is lowered 
passing p t h and p cr in succession. As for the lower section 
(7* < 7*) of L e = 1.8 diagram, upon increasing density, 
original homogeneous LL phase becomes destabilized and 
develops into the proposed ML phase where spatial het- 
erogeneity develops, until finally the RG phase takes the 
lead. At the crossing point, i.e. 7* =7*, RG melts into 
LL state without going via any intermediate phase. In 
the case of L e = 1.2, as we go across the phase boundary 
Pco by increasing density, the CO state would transform 
into the MS phase whereby the distinction between the 
two types of arrested states with different mechanisms of 
localization is recovered. 

The effective bond stiffness can be varied by manipu- 
lating the crosslinking and/or bundling properties or by 
changing the temperature. Some general features can be 
extracted from the presented state diagrams, that is high 
elastic stiffness tends to (1) stabilize the LL state and (2) 
facilitate the crossover to glassy behavior. The first ef- 
fect is quite explicit in L e = 1.8 case: as 7* increases, 
the elastic-nonlinearity-induccd ML phase evolves into 
LL state with a single stable au q when p < p cr , whereas 
purely RG state develops into the MS state when p > p cr . 
The second effect is clear in L e = 1.2 case; starting either 
from distinct LL state or from the MS phase, the system 
would end up with CO behavior as long as 7* transcends 
the phase boundary 7* (/c) (inversion of p co {l*))- 

Actually the parameter-modulated transformation of 
the phase behavior in terms of the order parameter a 
can be directly detected in a ta gged versus a n eighbor plots, 



which explicitly show the emergence and disappearance 
of, as well as transitions among, various fixed points un- 
der the parameter control. In other words, the phase 
boundaries essentially indicate switching between differ- 
ent fixed point structures of the self-consistent equations. 
For example, the p co boundary marks the disappearance 
of the lowest-a fixed point: in LL— s-CO case, modest dis- 
continuity in a value arises from stability shift to a newly 
established fixed point in proximity; while a considerable 
jump in a value is observed in MS— >CO case, since no 
intermediate fixed point develops, the high-a fixed point 
becomes the only stable attractor. Across the 7 t *^(p) (in- 
version of pth{l*)) boundary from below, the initially 
bifurcation-generating unstable fixed point becomes sta- 
bilized by enhanced stiffness. 

In addition to the information obtained from the me- 
chanical stability, the integrated thermodynamic charac- 
teristics are more informative of the glassy aspects. In the 
direction of vitrification, kinetic laboratory glass transi- 
tion is expected at pc in view of the landscape-dominated 
transport mechanism triggered at pa and the presence of 
extensively many possible frozen-in states^—. Whereas 
the Kauzmann density pk , at which the configurational 
entropy ceases to be extensive and the glassy configu- 
rations arc no longer metastable, indicates a thermody- 
namic transition that ultimately may underlie the kinetic 
arrest. Though hard to achieve on practical timescale, 
Pk does provide a mean-field estimate of how dense the 
liquid can be below which a glass transition would be 
forced to intervene to avoid the entropy crisis^. 

6. Bonded- fraction dependence of mechanical and 
thermodynamical properties 

So far, we have implicitly assumed that all the nearest- 
neighbor pairs arc bonded, i.e. the network is fully con- 
nected and there are no free beads at all. In biolog- 
ical fact, however, apart from the fibrous cytoskeletal 
network, there also exists a colloidal suspension of pro- 
tein molecules (including detached ABPs, recycled actin 
monomers, etc.) that contributes equally, if not more, to 
the crowding interior of a cell, and thus to the excluded 
volume effect. We mimic such a suspension of molecules 
simply as a collection of free beads, in which is immersed 
the nonlinear elastic fiber network which is anchored on 
the bonded beads. In our mean-field context, the fraction 
of bonded beads against the free ones is equivalent to the 
probability for a nearest-neighbor pair to be bonded. We 
assign an independent parameter Pj, £ (0, 1] to indicate 
this bonded fraction or network connectivity, and assume 
P, to be independent of the overall bead density and the 
effective bond stiffness to purify its influence. 

The self-consistent equation to determine a and the 
expressions of fn q and f g i are modified accordingly: 
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In equations for a and f g i, P b and (1 — P b ) lead the ef- 
fective potential between bonded and non-bonded pairs, 
respectively; and the influence due to bonded fraction 
change on the cell constraint term (i.e. the second order 
term in the effective potential expansion) is contained 
in the self-consistently determined a g i. As for fu q , the 
bonded fraction not only modifies the bonding correc- 
tion to the HS interaction, it also separates the bonded 
from non-bonded contributions to the entropy cost: for 
free particles InpA 3 — 1 should suffice to describe the 
density dependence of the entropy cost, whereas we use 
| In (aiiqA? /ite) — 1 for the bonded beads. Notice the fact 
that as bonds melt (i.e. P b decreases), increasing transla- 
tional symmetry would lower the entropy cost to localize 
the density waves. We shall show that the logarithmic 
dependence on au q used here could at least qualitatively 
incorporate this feature. Moreover, the bonding entropy 
due to various choices of bonded pairs among nearest 
neighbors is not explicitly included, since only the differ- 
ence between f g i and fu q matters for current purposes. 
It is easily seen that as P b — > 1 our earlier expressions for 
a fully connected network (Eqs. © (JT]) © ) are recovered. 

We show in Fig. [T7] and Fig. [15] the state diagrams 
at P b = 0.8 and P b = 0.5 for both typical networks. 
The corresponding contour maps of au q and a g i solu- 
tions (not shown here) indicate that the behavior of the 
liquid-like localization strength is not qualitatively af- 
fected by changing the bonded fraction, and the glassy 
solutions are almost quantitatively intact. Yet the tran- 
sition boundaries are significantly shifted as the bonded 
fraction varies. For the floppy network (Fig. IT7|) . as P b 
decreases, the slower increase of pth with 7* (which im- 
plies a lower destabilization density for the liquid-like so- 
lution) yields a shrinking MS region and enlarged RG 
and possibly ML regions. This behavior can be under- 
stood as arising from the fact that weaker network con- 
nectivity makes the stabilization via "bond trapping" less 
efficient. In the tense network (Fig. ITS)) , as the bonded 
fraction drops, both p t h and Pco boundaries shift upward 
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FIG. 17: The state diagrams of a typical floppy network (L e = 
1.8) with different bonded fraction. pK (red) and pa (green) 
are defined as former, (a) P b = 0.8; (b) Pj, = 0.5: transition 
densities are absent due to the limited fty range shown here. 



resulting in an extension of both MS and RG regions into 
a higher 7* regime. This observation implies a higher 
bond stiffness is needed to stabilize the liquid-like solu- 
tions so as to trigger the crossover to glassy behavior. As 
for the dynamical transition density P a (i.e. the critical 
density p cr analyzed in earlier sections), it modestly in- 
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FIG. 18: The state diagrams of a typical tense network [L e = 
1.2) with different bonded fraction, (a) P b = 0.8; (b) P b = 0.5. 
pG and pk become lower as the bonded fraction decreases. 



creases with lowering P], in the tense network (rises from 
0.91 to 0.97 as Pb drops from 1.0 to 0.5) while it remains 
constant (~ 1) in the floppy network; concomitantly a a 
halves its value (115 — > 55) in the tense network while it 
stays the same (~ 50) for the floppy case. 

We next examine the variation in thermodynamics due 
to the change of bonded fraction. In both networks, pc 
and pk are found to persist upon decrease in P\ (with 
the order of pa < Pg < Pk maintained). In the tense 
network, pq and pk develop moderate 7* dependence 
and shift toward lower density as Pb decreases, as shown 
in Fig. 1181 In the floppy network, pq and pk emerge 
at considerably higher 7* as Pb drops indicating greater 
difficulty in stabilizing LL motion, yet become insensitive 
to /?7-value thereafter, as seen in Fig. [T7Ta). When Pb is 
further lowered to 0.5 (i.e. network being half-connected) 
transition densities are absent due to the limited range 
of /?7 shown here (see Fig. [TTTb)) and would reappear if 
we extend /37 sufficiently. 

In contrast, if we use In pA 3 — 1 for both the bonded 
and non-bonded contributions to the entropy cost, a dra- 
matic change in pq and px is found (not shown here): 
when Pb = 0.8, in both networks, px is barely above 
Pa while pa is entirely skipped; if Pb is further lowered 
to 0.5, then fu q > f g i for all p > pa in MS region, in- 
dicating a negative configurational entropy which is not 



physically meaningful. Actually in our model the highly- 
localized glassy motion is insensitive to the degree of net- 
work connectivity since (/3V£/ del (a g i ; p > p A ;^*,L e )) ~ 
(PV$(a gl ;p > p A )) and (V 2 (3V^J del (a g r, p > 
p A ;i*,L e )) ~ (V 2 fiVg f J (a g r, p > pa)), thus such a 
significant drop in transition densities results from an 
enhanced attractive interaction in the liquid-like phase 
due to stronger thermal fluctuations (smaller au q ) in- 
duced by reduced bond constraints (lower Pb). This 
energetic enhancement in fu q is balanced, partly, by 
the decrease in entropic cost (to localize density wave) 
when I ln(a;j g A 2 /7re) — 1 is used, thereby mitigating the 
bonded-fraction modulation upon the transition densi- 
ties, and resulting in a persistent transition possibility 
over a large Pb range. 

In sum, the overall tendencies are clear: decreases in 
the bonded fraction lower the transition densities pa and 
Pk implying that the system becomes less capable of re- 
configuring (or easier to become glassy) upon weaken- 
ing of the network connectivity; in this sense, the model 
nonlincar-clastic-bonded interaction encourages liquid- 
like motion and thereby facilitates more efficient struc- 
tural rearrangements. 



IV. CONCLUSIONS 

In this paper, we have modeled the cytoskclcton as 
an amorphous network of rigidly cross-linked nonlinear- 
clastic bonds that become tense beyond an intrinsic onset 
length and buckle otherwise. We study the equilibrium 
mechanical properties of the model system within the 
framework established by the self-consistent phonon the- 
ory and the free energy functional formulation. 

We have obtained an initial understanding of the phys- 
ical behavior via the calculation of several representative 
thermodynamic quantities, and by examining the state 
diagram of typical systems. Diverse mechanical proper- 
ties of a generic cytoskclcton can be recognized by ana- 
lyzing the featured phases and possible transitions: the 
permanent network structure excludes a completely er- 
godic fluid phase, whereas the nonlinearity in the elas- 
tic interaction induces spatial heterogeneity that exhibits 
through a "martensitic-like" phase with domains of ori- 
ented distortions. The probable coexistence of the liquid- 
like and glassy behavior implies the capability of mak- 
ing structural rearrangements with varying agility in re- 
sponse to mechanical stimuli. The effective bond stiff- 
ness tends to stabilize the liquid-like state and facili- 
tates its crossover to glassy behavior, whereas the rel- 
ative position of the elasticity onset with respect to the 
nearest-neighbor shell dramatically modulates the transi- 
tion boundaries; the critical density may no longer mark 
a sharp transition in certain situations when crossover 
takes place. In sum, the elasticity onset length deter- 
mines all possible mechanical phases within a practical 
parameter range, while the bond stiffness decides which 
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transitions occur upon the variation of cross-link concen- 
tration (thus of the bead density). 

We further investigated how the bonded fraction or 
network connectivity modulates the phase boundaries as 
well as the thermodynamic terms of the transition densi- 
ties (pa and pk)- We found that decreasing the bonded 
fraction results in an upward shift of both p t h and p co 
boundaries, indicating the need for a higher bond stiff- 
ness to compensate for any loss in connectivity, so as to 
stabilize the liquid-like motion and to trigger its crossover 
to glassy behavior; on the other hand, the expanded 
multiple-solution phase region allows for a large stiffness 
range with extensive conhgurational degrees of freedom. 
As for the thermodynamics, the characteristic densities 
show little dependence on bond stiffness in floppy net- 
works; for tense networks, however, pa and Pk become 
lower upon enhanced bond stiffness, suggesting decreas- 
ing configurational degrees of freedom at a certain den- 
sity as bonds stiffen. Further, for a given density, possible 
glass transition takes place at a much lower stiffness in 
tense networks than in floppy ones, while for a certain 
bond stiffness, the model system becomes vitrified at a 
lower bead density as more bonds tense up. 

As exhibited clearly in the tense network, the loga- 
rithmic dependence of the entropy cost (for bonded in- 
teraction in the liquid-like phase) on particle localization 
strength adopted in our free energy functional contains 
the feature that as more bonds form, the kinetic glass 
transition would occur at a higher density, suggesting 
that the nonlinear-elastic-bonded interaction might help 
resolve local/steric constraints and facilitate escape from 
topological trapping, resulting in a more dense packing 



when stuck finally. Fleshing out this conjecture may re- 
quire coming to term with finite range consideration; we 
need to go beyond mean-field level and consider activated 
events among various metastable states (probably via 
"droplet relaxation" in a mosaic structure^) . In a biolog- 
ical sense, cells may prefer an interconnected structural 
skeleton, not only to maintain their architecture, but 
to realize more efficient structural rearrangements when 
necessary. As argued by Ingber, the tensed/prestressed 
hierarchical networks play a central role in producing a 
well-orchestrated multiscale mechanical response^. 

This work provides a general scheme to study macro- 
scopic mechanical phases in terms of the stability to lo- 
cal mechanical environment, and our current equilibrium 
model sets up a test field for further incorporated features 
for a more realistic model, in particular, the motorization 
effect that makes the system active, far from equilibrium 
and makes it physically distinct from an ordinary poly- 
mer network. It will be interesting to study the interplay 
of bond constraints and force-environment-sensitive mo- 
tors in maintaining the cell's architecture and modulating 
the transition behavior. We also plan to investigate the 
effect of spatial heterogeneity and visualize the structural 
rearrangements by combining analytic schemes with sim- 
ulation techniques. 
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